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Abstract 

Four  models  of  convective  and  radiative  heat  transfer  inside  tubular  solid  oxide  fuel  cells  are  presented  in  this  paper,  all  of  them  applicable  to 
multidimensional  simulations.  The  work  is  aimed  at  assessing  if  it  is  necessary  to  use  a  very  detailed  and  complicated  model  to  simulate  heat 
transfer  inside  this  kind  of  device  and,  for  those  cases  when  simple  models  can  be  used,  the  errors  are  estimated  and  compared  to  those  of  the  more 
complex  models. 

For  the  convective  heat  transfer,  two  models  are  presented.  One  of  them  accounts  for  the  variation  of  film  coefficient  as  a  function  of  local 
temperature  and  composition.  This  model  gives  a  local  value  for  the  heat  transfer  coefficients  and  establishes  the  thermal  entry  length.  The  second 
model  employs  an  average  value  of  the  transfer  coefficient,  which  is  applied  to  the  whole  length  of  the  duct  being  studied.  It  is  concluded  that, 
unless  there  is  a  need  to  calculate  local  temperatures,  a  simple  model  can  be  used  to  evaluate  the  global  performance  of  the  cell  with  satisfactory 
accuracy. 

For  the  radiation  heat  transfer,  two  models  are  presented  again.  One  of  them  considers  radial  radiation  exclusively  and,  thus,  radiative  exchange 
between  adjacent  cells  is  neglected.  On  the  other  hand,  the  second  model  accounts  for  radiation  in  all  directions  but  increases  substantially  the 
complexity  of  the  problem.  For  this  case,  it  is  concluded  that  deviations  between  both  models  are  higher  than  for  convection.  Actually,  using  a 
simple  model  can  lead  to  a  not  negligible  underestimation  of  the  temperature  of  the  cell. 

©  2007  Elsevier  B.V.  All  rights  reserved. 
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1.  Introduction 

SOFCs  are  devices  operating  at  temperatures  ranging  from 
800  to  1050  °C  for  state  of  the  art  materials.  Below  this  range, 
voltage  losses  due  to  ionic/electronic  resistivity  of  materials 
increase  noticeably  as  conductivity  grows  exponentially  with 
temperature  [1,2],  On  the  other  hand,  SOFCs  cannot  be  oper¬ 
ated  continuously  at  a  very  high  temperature,  say  1 100  °C,  as  this 
would  lead  to  a  considerable  decrease  in  performance,  probably 
caused  by  a  thermal  expansion  mismatch  between  electrodes 
and  electrolyte  [3].  Therefore,  the  management  of  heat  trans¬ 
fer  inside  a  solid  oxide  fuel  cell,  either  with  tubular  or  planar 
technology,  is  essential  in  order  to  guarantee  the  reliability  and 
long  life  demanded  by  the  market  to  this  sort  of  power  genera- 
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tion  devices.  Fig.  1  shows  the  amount  of  energy  released  and/or 


consumed  inside  an  SOFC  fed  with  natural  gas  as 
operating  voltage  and  for  different  pressures. 

;  a  function  of 

Three  reactions  are  considered  to  take  place:  hydrogen  oxida¬ 
tion,  Eq.  (1),  methane  reforming,  Eq.  (2),  and  carbon  monoxide 

shifting,  Eq.  (3). 

h2  +  |o2  -*  H20 

(1) 

CH4  +  H2O  -*  3H2  +  CO 

(2) 

CO  +  h2o  -+  h2  +  co2 

(3) 

The  net  amount  of  heat  released  according  to  Fig.  1  must  be 
evacuated  from  inside  the  cell  by  the  air  mass  flow,  which  is 
supplied  well  in  excess  with  respect  to  the  stoichiometry  of  Eq. 
(1).  Thus,  under  normal  operating  conditions,  only  15-20%  of 
the  air  is  used  to  oxidize  the  fuel. 

This  work  deals  with  heat  transfer  characterization  and  mod¬ 
elling  inside  tubular  SOFCs,  particularly  applied  to  a  1.5  m  long 
Siemens  Westinghouse  cell  with  100  W  rated  power  for  ambi- 
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Nomenclature 

A  cross-sectional  area  [m2] 

ASP  air  supply  pipe 

CT  matrix  of  heat  transfer  coefficients 

D  diameter  [m] 

/  friction  factor 

F  view  factor 

G  vector  of  generation  terms 

hcv  convective  heat  transfer  coefficient  [W  m-2  K- 1  ] 

J  radiosity  [W  m-2] 

k  thermal  conductivity  [W  m-1  K-1] 

L  cell  length  [m] 

m  mass  flow  [kg  s_  1  ] 

n  slice  number 

Nu  Nusselt  number 

Pr  Prandtl  number 

q  heat  flow  per  unit  area  [W  m-2] 

Q  total  heat  flow  [W] 

Re  Reynolds  number 

S  wall  surface  [m2] 

T  temperature  [K] 

T  vector  of  temperatures  [K] 

Uf  fuel  utilization  factor  [%] 

x  molar  fraction 

Gentry  thermal  entry  length  [m] 

Greek  symbols 
a  fluid  property 

e  emissivity 

lx  dynamic  viscosity  [p.P] 

er  Steffan-Boltzmann’s  constant 

Subscripts 
an  anode 

ca  cathode 

h  hydraulic 

J  relative  to  radiosity 

lam  laminar 

t  iteration  step 

tur  turbulent 

T  relative  to  temperature 

wall  property  of  a  wall/surface 

Superscripts 

*  normalized  view  factor 


ent  pressure  operation.  Fig.  2.  More  precise  geometric  data  of 
this  technology  can  be  found  in  reference  [4]  and  is  shown  in 
Table  1. 

As  said  before,  heat  released  in  Eq.  (1)  is  evacuated  from 
the  electrodes/electrolyte  solid  structure,  also  known  as  PEN 
from  Positive  Electrolyte  Negative,  mainly  by  convection  but, 
in  the  case  of  the  tubular  technology  shown  in  Fig.  2,  radiation 
between  PEN  and  air  supply  pipe  also  plays  an  important  role. 
Fig.  3  shows  the  proportion  of  total  heat  transfer  which  takes 


Fig.  2.  Reference  geometry. 


place  by  convection  and  radiation.  It  can  be  concluded  that  each 
one  of  these  heat  transfer  mechanisms  is  dominant  at  a  different 
part  of  the  cell:  convection  for  the  first  third  of  it  and  radiation 
from  that  point  to  the  exhaust  section.  Although  this  distribution 
is  not  constant  and  may  vary  according  to  operating  conditions, 
it  is  clear  that  these  two  phenomena  must  be  very  well  described 
when  developing  a  model  of  performance  suitable  for  tubular 
SOFCs. 

The  work  is  divided  in  two  parts.  First,  a  model  for  describing 
convective  heat  transfer,  based  on  a  local  evaluation  of  trans¬ 
fer  coefficients,  is  proposed.  This  model  is  later  simplified  and 
the  loss  of  accuracy  evaluated.  Secondly,  two  models  of  radia- 


Table  1 

Reference  geometry 


Length  [m]  1.5 

Anode  outer  diameter  [mm]  22 

Anode  thickness  [pm]  100 

Electrolyte  thickness  [tun]  40 

Cathode  thickness  [mm]  2.2 

Metallic  interconnection  thickness  [p.m]  85 
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Fig.  3.  Heat  transfer  by  convection/radiation. 


five  heat  transfer  are  proposed.  One  of  them  is  a  simple  model, 
extensively  used  in  previous  works,  for  radiative  exchange  in  the 
radial  direction  which  is  based  on  the  hypothesis  of  infinite  walls. 
Then,  a  complete  model  of  radiation  is  described.  This  second 
model  considers  radiative  exchange  in  all  directions,  radially  and 
obliquely,  and  introduces  additional  complexity  to  heat  balance 
equations. 

Fig.  4  shows  the  discretization  of  the  cell  which  is  used  for  the 
heat  transfer  models.  The  cell  is  divided  axially  into  a  number  of 
slices  which  are  again  divided  radially  into  five  annular  volumes, 
cylindrical  for  the  inner  one,  called  elements. 

It  is  commonly  agreed  that  multidimensional  models  of  per¬ 
formance  of  an  SOFC  are  based  on  a  decoupled  solving  strategy 
for  temperatures  and  composition.  In  other  words,  an  iterative 
method  is  used  that  solves  thermal  and  electrochemical  models 
consecutively  [5].  Firstly,  temperatures  are  calculated  by  solving 
the  system  of  linear  equations  resulting  from  local  heat  balance 
equations.  Compositions  and  current  density  remain  constant  at 
all  the  slices  and  elements.  This  system  of  equations  is  shown 
in  Eq.  (4)  where  T  stands  for  local  temperatures,  CT  for  heat 


i-1  SLICE  i  i+1 


transfer  coefficients  and  G  for  generation  terms: 

CT  T  =  G  (4) 

This  temperature  field  is  then  used  as  an  input  to  solve  the 
electrochemical  model  and  obtain  compositions,  current  density, 
voltage  losses,  reaction  rates,  etc. 

As  stated  above,  this  work  is  aimed  at  describing  radiative  and 
convective  heat  transfer  equations  involved  in  heat  balance  local 
equations.  In  other  words,  models  will  be  presented  to  calculate 
coefficients  included  in  CT,  Eq.  (4). 

2.  Convective  heat  transfer:  model  description 

Convective  heat  transfer  is  described  by  Newton’s  law  of 
cooling: 

<7conv  =  ^cv(7wall  —  Tgas)  (5) 

where  hcv  is  the  convective  heat  transfer  or  film  coefficient.  Cal¬ 
culating  this  coefficient  accurately  is  the  key  task  to  obtaining 
a  precise  heat  transfer  model.  The  following  lines  describe  a 
step-by-step  procedure  to  obtain  hcv: 

1.  Calculation  of  fluid  properties:  viscosity  and  thermal  con¬ 
ductivity. 

2.  Calculation  of  Reynolds  number  from  fluid  properties  and 
duct  geometry. 

3.  Calculation  of  flow  regime  from  Reynolds  number. 

4.  Calculation  of  Nusselt  number  and,  consequently,  convective 
heat  transfer  coefficient. 

2.1.  Properties  of  gases 

Two  gas  properties  are  needed  to  evaluate  hcv  at  each  wall: 
dynamic  viscosity  and  thermal  conductivity.  However,  two 
major  difficulties  emerge  that  make  it  complex.  Firstly,  the  lack 
of  experimental  data  at  the  very  high  temperature  and  particu¬ 
lar  composition  which  are  characteristic  of  SOFCs,  especially 
for  conductivity.  Besides,  these  two  properties  do  not  depend 
on  composition  proportionally  because  of  the  very  different 
behaviours  of  most  of  the  seven  species  considered  in  the  mix¬ 
ture:  methane,  carbon  monoxide,  carbon  dioxide,  hydrogen, 
water  vapour,  oxygen  and  nitrogen.  In  other  words,  they  cannot 
be  calculated  with  a  general  expression  like: 

«gas(T)  =  •  X{  (6) 

which  is  valid  for  specific  heat  for  example,  a,-  and  x,-  stand  for 
property  and  molar  fraction  of  pure  components,  respectively, 
in  Eq.  (6). 

The  first  step  is  calculating  the  dynamic  viscosity  of  pure 
components  as  a  function  of  temperature.  For  this  purpose,  a 
fifth  order  polynomial  given  in  reference  [6]  is  used,  Eq.  (7). 
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Table  2 

Coefficients  for  dynamic  viscosity  and  conductivity  calculations 


CH4  CO  C02  h2  h2o  o2  n2 


aHA 

an5 


ak,l 

Ok,2 


ait, 5 
Ok,6 


— 9.9989  -4.9137 

529.37  793.65 

-543.82  875.90 

548.11  883.75 

-367.06  -572.14 

140.48  208.42 

-22.920  -32.298 

0.4796  -0.2815 

1.8732  13.999 

37.413  -23.186 

-47.440  36.018 

38.251  -30.818 

-17.283  13.379 

3.2774  -2.3224 


-20.434  15.553 

680.07  299.78 

-432.49  -244.34 

244.22  249.41 

-85.929  -167.51 

14.450  62.966 

-0.4564  -9.9892 

2.8888  1.5030 

-27.018  62.892 

129.65  -47.190 

-233.29  47.763 

216.83  -31.939 

-101.12  11.972 

18.698  -1.8954 


-6.7541  -1.6918 

244.93  889.75 

419.50  -892.79 

-522.38  905.98 

348.12  -598.36 

-126.96  221.64 

19.591  -34.754 

2.0103  -0.1857 

-7.9139  11.118 

35.922  -7.3734 

-41.390  6.7130 

35.993  -4.1797 

-18.974  1.4190 

4.1531  -0.2278 


1.2719 

771.45 

-809.20 

832.47 

-553.93 

206.15 

-32.430 

-0.3216 

14.810 

-25.473 

38.837 

-32.133 

13.493 

-2.2741 


Coefficients  can  be  found  in  Table  2. 


A  similar  expression  is  used  for  thermal  conductivity: 


3.  Anodic  duct.  As  shown  in  Fig.  5,  the  geometry  of  the  anodic 
duct  is  rather  complex  and  its  hydraulic  diameter  is  calculated 
from  Eq.  (12): 

Cross-sectional  area  4  —  jt 

Oh  =  4—— - — -  =j  - Oan out-  (12) 

Wet  perimeter  n 


Next  step  after  calculating  pure  components  properties  is  to 
evaluate  those  of  the  mixture.  According  to  the  work  by  Todd  and 
Young  [6],  where  several  methods  used  to  calculate  fluid  prop¬ 
erties  are  studied  and  compared,  the  most  accurate  methods  are 
Reichenberg’s  and  Wassiljevas’s  for  viscosity  and  conductivity, 
respectively.  Expressions  for  these  methods  are  rather  complex 
and  will  not  be  quoted  here;  however,  readers  wishing  to  know 
full  expressions  are  referred  to  [6]. 

2.2.  Reynolds  number 

Reynolds  number  evaluates  the  ratio  of  viscosity  to  inertia  of 
a  gas  or  liquid  stream  and  is  representative  of  the  flow  regime  of 
a  stream.  The  following  general  expression  of  Re  inside  a  duct 
is  used: 


2.3.  Flow  regime  and  entry  length 

Three  cases  are  considered  when  studying  the  flow  regime  of 
each  stream  according  to  Re: 

1 .  Laminar  flow,  thermal  entry  length. 

2.  Laminar  flow,  fully  developed. 

3.  Turbulent  flow,  fully  developed. 

Although  many  authors  claim  that  fully  turbulent  flow  is  not 
guaranteed  below  10,000,  a  critical  value  of  Re  is  set  at  2300 
according  to  [1,7,8],  For  laminar  flow,  the  following  thermal 
entry  length  is  considered  [7] 

Xentry, lam  =  0.0575  Re  Pr  Dh  (13) 

where  thermal  and  hydrodynamic  boundary  layers  develop 
simultaneously  due  to  a  Prandtl  number  close  to  unity,  Pr  =  0.7. 
For  turbulent  flow,  the  evaluation  of  entry  length  is  more  com- 


(9) 


where  A  is  the  cross-sectional  area,  m  the  mass  flow  and  Dh  is 
the  hydraulic  diameter  or  characteristic  length  of  the  duct.  This 
diameter  is  different  for  each  duct: 


1.  Air  supply  pipe.  Hydraulic  and  geometric  diameters  are  the 
same. 

Dh  =  Daspjti  (10) 

2.  Annular  duct.  The  hydraulic  diameter  is: 

(ID 


Dh  =  Dca  in  -  DASP.out 
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Table  3 

Nusselt  number 


x*  =  x/Re  Pr  -D  Nux 


0.001  16.8 

0.002  12.6 

0.004  9.6 

0.006  8.25 

0.010  6.8 

0.020  5.3 

0.050  4.2 


oo  3657 

Laminar  entry  length. 

plex  as  the  transition  from  undeveloped  to  fully  developed  flow 
is  thought  to  take  place  somewhere  inside  the  range: 

1Q  ^  Xentryturb  <  g()  (14) 

Dh 

The  ratio  of  total  length  to  hydraulic  diameter  of  the  ducts 
involved  in  this  study  ranges  from  150,  air  supply  pipe,  to  280, 
annular  duct.  For  these  geometric  “boundary”  conditions,  the 
turbulent  entry  length  is  considered  negligible. 


Fig.  6.  Correlations  for  turbulent  flow  [10]. 


2.4.  Nusselt  number:  convective  heat  transfer  coefficient 


The  heat  transfer  model  presented  in  this  work  is  the  evolu¬ 
tion  of  a  multidimensional  SOFC  model  previously  developed 
by  the  authors  [5,9]  which  predicts  temperatures  and  composi¬ 
tions  locally,  at  each  point  of  the  cell.  Thus,  the  evaluation  of 
heat  transfer  coefficients,  i.e.  Nu,  must  include  a  dependence  on 
position,  and  consequently  composition,  inside  the  cell. 

Nusselt  number  is  calculated  from  correlations  fitted  to 
empirical  data,  either  by  means  of  mathematical  expressions 
or  tables.  Most  of  these  expressions  give  values  for  average  Nu 
along  the  whole  duct  and  only  a  few  of  them  are  applicable  to 
local  studies.  The  latter  will  be  used  in  this  work. 

Values  of  Nu  in  Table  3  [8]  are  used  for  the  laminar  entry 
length.  Interpolation  is  used  between  given  non-dimensional 
positions  x* .  The  last  value  of  Nu  in  Table  3  corresponds  to 
the  value  for  fully  developed  laminar  flow  which  holds  constant 
at  3.657  for  any  value  of  Re  or  x. 

For  Re>Rec,  flow  is  turbulent  and  Gnielinski’s  equation  is 
used  to  evaluate  Nu,  Eqs.  (15)  and  (16)  [7],  where/stands  for 
friction  factor  in  Eq.  (16). 


(f/8)(Re  -  1 000]  Pr 
'  1  +  12.7 V^P)(Er2/3  -  1) 


■*(?) 


0.791n(/?e)  —  1.64 


(15) 

(16) 


Gienelinski’s  equation  is  applicable  to  Re  >  2300, 
0.5  <  TV <2000  and  L>Dh.  The  complexity  of  Eqs.  (15) 
and  (16)  leads  some  authors  to  employ  simpler  equations  like 
Colburn’s,  shown  in  Eq.  (17),  which  is  valid  for  Re>  10,000, 
0.7  <Pr<  160  and  L  >  1  ()l\ .  This  correlation  is  easier  to 
evaluate  but  can  lead  to  errors  as  high  as  20%  as  shown  in  Fig.  6 
[10].  In  addition,  the  uncertainty  in  the  range  of  Re  from  2300, 


end  of  laminar,  to  10,000,  Eq.  (17)  validity,  is  significantly 
higher  than  when  using  Gnielinski’s  correlation. 

Nu  =  0.023  Re4/5Pr1/3  (17) 


Finally,  the  convective  heat  transfer  coefficient  is  evaluated 
directly  from  the  value  of  Nu  through  the  following  equation: 


/lCv 


Nuk 


(18) 


The  method  described  above  is  based  on  substituting  Dh  by 
x,  axial  position,  in  Eq.  (15)  and  gives  a  value  of  Nu  which 
is  not  the  real  value  of  local  Nu.  Fig.  7  shows  that  there  is  an 
underestimation  with  respect  to  the  actual  local  Nu.  However, 
the  error  made  when  using  this  approach  is  small  and  does  not 
increase  substantially  the  uncertainty  in  the  value  of  heat  transfer 
coefficients. 


3.  Radiative  heat  transfer:  model  description 

Previous  works  by  the  authors  have  shown  that  radiation 
involves  not  only  heat  exchanged  between  solid  walls  but 
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between  a  solid  wall  and  certain  gases.  However,  the  latter 
is  out  of  the  scope  of  this  work  as  it  is  only  relevant  under 
abnormal  operation  inside  the  cell,  i.e.  high  current  density 
[5,9]. 

Two  models  are  to  be  presented  in  the  following  lines.  One 
of  them  is  a  so-called  radial  model  and  is  based  on  the  hypoth¬ 
esis  of  infinite  coaxial  cylinders.  It  will  be  called  the  simple 
model.  A  more  complex  model  will  be  presented  later,  which 
accounts  for  oblique  radiation  between  axially  adjacent  cells. 
This  model,  called  complex  model,  makes  it  necessary  to  intro¬ 
duce  a  new  unknown  for  each  wall  involved  in  the  radiative 
exchange. 

In  both  models,  exchanging  walls  are  considered  to  be  grey 
walls  satisfying  Kirchhoff’s  laws. 


3.1.  Simple  model:  infinite  parallel  exchanging  walls 

This  simple  model  assumes  exchanging  walls  to  be  infinitely 
long  and,  thus,  only  radiation  between  walls  of  the  same  slice. 
Fig.  4,  is  considered.  In  addition,  this  hypothesis  reduces  the 
radiative  exchange  to  the  walls  of  the  annular  cathodic  duct  as 
long  as  wall  temperatures  are  considered  to  be  constant  at  each 
slice.  In  other  words,  for  the  inner  duct  of  the  air  supply  pipe  and 
the  anodic  duct,  the  net  amount  of  heat  exchanged  by  radiation 
between  walls  is  zero. 

Eq.  (19)  defines  the  net  amount  of  heat  exchanged  by  radia¬ 
tion  between  the  outer  wall  of  the  supply  pipe  ASP  and  the  inner 
wall  of  the  cathode. 


'  (^cathode  -  ^ASp) 


(1/  ^cathode)  +  (Ocathode/^ASP)((l /£ASp)  —  1) 

X  7T  •  -Dcathode  '  Ax 


The  view  factor  between  inner  and  outer  surfaces  of  the  duct 
is  1  for  Fas  Pxat  and  the  ratio  of  diameters  for  FcatiASP- 

The  radiative  flow  calculated  from  Eq.  (19)  forms  part  of  the 
linear  system  of  equations  used  to  solve  the  thermal  model  of  the 
fuel  cell,  Eq.  (4),  where  T  is  a  vector  containing  temperatures 
to  be  calculated,  five  for  each  slice  of  the  cell.  CT  is  a  matrix 
containing  temperature  coefficients  for  each  heat  balance  equa¬ 
tion  and  G  stores  known  heat  fluxes  like  heats  of  reaction.  This 
system  of  equations  has  5 N  unknowns,  N  being  the  number  of 
slices  that  the  cell  is  divided  into. 

Finally,  Eq.  (19)  must  be  linearized  before  added  to  Eq.  (4). 
To  do  so,  the  following  approach  is  used: 


T)  =  4  -  7j_i  ■  T,  =  4  -  +  ^-  y  .  Tl  (20) 

where  t—  1  and  t  are  the  previous  and  current  steps  of  the 
iterative  process  and  Tm  is  the  arithmetic  mean  temperature 
of  cathode  and  air  supply  pipe.  This  approach  is  accurate 
enough  for  temperature  differences  between  walls  below  50  °C 
[11]. 


3.2.  Complex  model:  oblique  radiation 

The  previous  model  does  not  consider  oblique  radiation 
between  adjacent  slices  so  the  radiative  flow  leaving  the  surface 
through  both  annular  ends  of  each  slice  is  dismissed.  As  dis¬ 
cussed  later,  this  is  only  acceptable  for  a  low  number  of  slices, 
i.e.  high  slice  length.  When  this  is  not  the  case,  the  amount 
of  radiation  not  considered  for  heat  transfer  balance  equations 
can  be  not  negligible  and  lead  to  high  deviation  in  wall  tempera¬ 
ture.  In  such  a  situation,  a  more  complex  description  of  radiative 
exchange  is  needed. 

Unfortunately,  the  complexity  of  the  model  which  includes 
oblique  radiation  is  significantly  higher  than  that  of  the  simple 
model.  The  radiative  exchange  inside  a  multiple  wall  enclosure 
cannot  be  evaluated  from  temperatures  exclusively,  and  a  new 
set  of  parameters  is  needed.  In  fact,  each  wall  of  the  domain  is 
characterized  not  only  by  its  temperature  but  by  its  radiosity  as 
well,  the  latter  being  the  amount  of  thermal  radiation  leaving  the 
wall. 

At  a  certain  enclosure,  the  heat  flow  leaving  a  surface  is 
determined  by  the  following  equation: 


Grad.,  =  -  YsWjj  Si  (21) 

where  Jj  is  the  radiosity  of  surface  i,  Fq  the  view  factor  between 
walls  i  and  j,  respectively,  and  Sj  is  the  area  of  wall  i.  In  order  to 
evaluate  radiosities,  an  equation  for  each  surface  is  needed,  Eq. 
(22),  which  links  temperature  and  radiosity  at  each  wall  of  the 
enclosure: 


Ji  =  SioT*  +  ^T(l  -  sfiFijJj  (22) 

V; 

where  e,-  is  the  surface  emissivity  and  a  is  the  Stefan-Boltzman 
constant.  Again,  and  in  order  to  linearize  the  system  of  equations 
to  be  solved,  the  power  of  temperature  to  the  fourth  is  expressed 
in  the  following  manner: 

Ti=Tlt-\TUt  (23) 


where  i  and  t  stand  for  surface  and  iteration  step,  respectively. 

Eqs.  (21)  and  (22)  are  introduced  in  the  system  of  equations 
described  in  Eq.  (4)  to  form  a  new  one: 


Tctt  CTjl  M  _  Tgt1 
[cjj  CJ]\  [/J  ~~  [gjJ 


(24) 


Clr  stores  the  coefficients  related  to  temperatures  for  the  heat 
balance  equations  at  each  domain  of  each  slice.  C7j  stores  the 
coefficients  related  to  radiosity  for  the  same  heat  balance  equa¬ 
tions;  in  this  case,  the  coefficients  are  included  in  the  calculation 
of  radiative  heat  flows  according  to  Eq.  (21).  C.Jj  and  CJ\  store 
the  coefficients  related  to  temperature  and  radiosity,  respectively, 
in  Eq.  (22),  radiosity  balance  at  each  wall.  Independent  terms 
Gx  and  Gj  have  been  reorganized  conveniently. 
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Fig.  8.  Local  Re  inside  ducts  for  V=0.45  V. 


4.  Convective  heat  transfer:  model  results 

As  said  in  the  introductory  section  of  this  work,  convective 
heat  transfer  coefficients  cannot  be  taken  as  constant  along  the 
cell  tube.  In  fact,  as  shown  in  Fig.  8,  Re  varies  significantly 
from  the  entrance  to  the  exhaust  section  of  any  of  the  three 
ducts  inside  the  stack,  especially  at  the  anode  as  a  consequence 
of  the  rapid  increase  in  temperature  and  mass  flow.  Convective 
heat  transfer  coefficients  are  depicted  in  Fig.  9.  It  can  be  seen 
that,  for  the  operating  conditions  considered,  0.45  V  at  80%  fuel 
utilization,  laminar  flow  is  only  found  at  the  air  supply  pipe. 
In  this  case,  it  is  also  easy  to  identify  the  thermal  entry  length, 
where  convective  transfer  coefficients  decrease.  Turbulent  flow 
is  found  at  the  other  two  ducts.  In  addition,  at  the  anodic  duct, 
it  must  be  noticed  that  h  varies  significantly  due  to  two  main 
causes:  first,  a  strong  temperature  variation  along  the  first  twenty 
centimetres  of  the  tube;  second,  an  increasing  mass  flow  from 
the  entrance  to  the  exhaust. 

However,  the  values  of  h  found  at  Fig.  9  are  not  usual  for  nor¬ 
mal  operating  fuel  cells,  as  practical  voltages  are  around  0.62  V. 
For  these  conditions,  flow  is  laminar  at  all  three  ducts,  which 
show  thermal  entry  lengths.  Heat  transfer  coefficients  for  this 
case,  shown  in  Fig.  10,  are  closer  to  reality. 


At  this  point,  a  comparison  between  both  assumptions,  con¬ 
stant  or  variable  heat  transfer  coefficients,  is  mandatory  in  order 
to  assess  its  impact  on  fuel  cell  modelling.  First,  the  model  has 
been  run  for  0.65  and  0.35  V,  80%  fuel  utilization,  900  °C  oper¬ 
ating  temperature  and  3  bar  pressure,  keeping  intake  fuel  and 
air  flows  equal  for  both  cases.  Results  are  shown  in  Table  4. 
The  first  conclusion  drawn  is  that,  from  a  performance  point  of 
view,  there  is  no  need  to  consider  variations  of  transfer  coeffi¬ 
cients  inside  the  cell  as  its  influence  on  power,  current  density 
and  other  relevant  parameters  of  the  cell  is  hardly  one  percent 
of  their  total  value.  Only  mean  temperatures  are  affected  at  low 
current  density  due  to  the  effect  of  thermal  entry  lengths  being 
considered  in  the  detailed  model.  This  effect  has  been  further 
studied  by  running  the  simplified  model  with  a  different  set  of 
initial  and  boundary  conditions,  rightmost  column  in  Table  4. 
In  this  case,  the  air  flow  entering  the  cell  has  been  reduced  in 
order  to  keep  the  mean  temperature  around  900  °C,  which  is 
the  temperature  predicted  by  the  detailed  model,  and  80%  fuel 
utilization.  Thus,  if  the  mean  temperature  were  to  be  kept  con¬ 
stant,  a  little  increase  in  air  utilization  from  9  to  10%  for  the  base 
case  would  arise.  However,  this  difference  tends  to  decrease  with 
increasing  current  densities  as  flow  becomes  turbulent  in  every 
duct. 

An  internal  study  has  also  been  carried  on  and,  for  the  cases 
shown  in  Table  4,  a  comparison  between  internal  temperature 
fields  has  been  done.  Results  at  0.65  V  operating  voltage  are 
shown  in  Fig.  1 1  for  all  three  cases  mentioned  above.  It  can 
be  seen  that  predicted  temperatures  are  significantly  different  at 
both  ends  of  the  cell  and  that  the  general  shape  of  the  curves  are 
similar  elsewhere.  This  behaviour  was  expected  and  is  caused 
by  the  miscalculation  of  heat  transfer  coefficients  at  the  thermal 
entry  lengths. 

Fig.  1 1  shows  that  a  difference  of  around  30  °C  exists  between 
all  three  models  at  different  parts  of  the  cell  and,  however, 
according  to  results  in  Table  4,  this  is  not  noticeable  when  ana¬ 
lyzing  the  global  performance  of  the  stack.  It  has  been  said 
previously  that  local  temperature  affects  the  electrochemical  per¬ 
formance  of  the  cell  and,  thus,  it  has  a  major  impact  on  current 
density.  This  effect  is  depicted  in  Fig.  12,  where  it  can  be  seen 
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Table  4 

Comparison  at  constant  fuel  and  air  intake  flows 

Variable  Nu  Constant  Nu 

0.65  V  0.35  V  0.65  V  0.35  V  0.65  V 

Current  density  [Am“2]  1745  5382  1740  5392  1758 

Power  [W]  94.6  157.4  94.3  157.4  95.3 

Fuel  utilization  [%]  79.8  80.3  79.5  80.5  79.7 

Temperature  [°C]  902(922)  898(928)  886(906)  900(928)  902(926) 

Maximum  temperature  shown  in  parentheses  in  the  mean  temperature  box. 


that  the  local  current  density  depends  strongly  on  temperature 
but,  as  long  as  the  mean  temperature  is  similar  for  all  three  cases, 
the  mean  current  density  is  equally  similar.  For  this  reason,  the 
global  performance  of  the  cell  is  not  affected. 

5.  Radiative  heat  transfer:  model  results 

Results  obtained  when  applying  the  simple  model  have 
been  published  in  previous  works  by  the  authors  [9]  and  other 
researchers  [4]  and  will  not  be  repeated  here.  However,  it  must 
be  said  that  they  have  been  considered  satisfactory  as  the  impact 


CATHODE 


AIR  SUPPLY  PIPE 


Fig.  13.  Radiative  exchange  of  cathode  inner  face  (wall  A). 

over  global  performance  of  using  the  simple  or  complex  models 
is  not  dramatic.  This  will  be  shown  later. 

Fig.  13  shows  the  complexity  of  the  radiative  heat  transfer 
when  applying  the  complex  model.  A  is  the  inner  wall  of  the 
cathode  of  slice  n  and  is  exchanging  heat  by  radiation  with  itself, 
the  inner  wall  of  the  cathode  of  adjacent  slices  n  +  1  and  n—  I 
(wall  C)  and  the  outer  wall  of  the  air  supply  pipe  of  the  same 
slice  (wall  B)  and  adjacent  slices  n+  1  and  n  —  1  (wall  D).  Radia¬ 
tive  exchange  with  slices  further  than  n  +  2  or  n  —  2,  depicted  in 
Fig.  13  with  dashed  lines,  is  considered  to  take  place  as  a  whole 
through  wall  E. 

Table  5  shows  the  calculated  view  factors  when  considering 
100  slices.  For  this  case,  if  the  radiation  exchanged  between 
A  and  walls  further  than  n  +  2  is  not  considered,  about  20%  of 
the  total  amount  of  radiation  emitted  by  A  is  neglected  in  the 
heat  balance.  In  addition,  this  percentage  of  radiation  lost  tends 
to  increase  when  the  number  of  slices  is  high  and  slices  get 
narrower.  Two  strategies  are  proposed  in  order  to  take  this  issue 
into  account. 

The  first  method  is  based  on  additional  exchange  equations. 
Radiation  between  slices  which  are  further  than  n  +  2  or  n  —  2 
can  be  added  until  the  equivalent  Fae  view  factor  is  smaller  than 
a  maximum  value  previously  established.  This  method  increases 
the  computation  time  and,  as  a  main  drawback,  it  can  generate 
a  badly  conditioned  exchange  matrix  C7j  or  CJj  where  some 
elements  are  too  small  compared  with  elements  on  the  main 
diagonal. 

The  second  approach  to  the  problem  is  based  on  neglecting 
the  radiation  lost  through  the  virtual  wall  E  and  normalizing  the 

Table  5 

View  factors  for  a  hundred  slices 

Fab  Faa  Fac  Fad  Fae 

0.4163  0.2095  0.0456  0.0319  0.1096 


Fig.  12.  Current  density  distribution  along  the  cell. 
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Table  6 

Comparison  at  constant  fuel  and  air  intake  flows 


Simple  model 

Complex  model 

0.65  V 

0.35  V 

0.65  V 

0.35  V 

Current  density 
[Am"2] 

1721 

5259 

1745 

5382 

Power  [W] 

93.3 

153.5 

94.6 

157.4 

Fuel  utilization  [%] 

78.7 

78 

79.8 

80.3 

Temperature  [°C] 

851(884) 

881(927) 

902(922) 

898(928) 

Maximum  temperature  shown  in  parentheses  in  the  mean  temperature  box. 


Table  7 

Comparison  at  80%  fuel  utilization  and  900  °C 


Simple  model 

Complex  model 

0.65  V 

0.35  V 

0.65  V 

0.35  V 

Current  density  [Am-2] 

1721 

5259 

1745 

5382 

Power  [W] 

93.3 

153.5 

94.6 

157.4 

Fuel  utilization  [%] 

78.7 

78 

79.8 

80.3 

Air  utilization  [%] 

14 

7.5 

9 

6.5 

Temperature  [°C] 

900(943) 

902(954) 

902(922) 

898(928) 

Maximum  temperature  shown  in  parentheses  in  the  mean  temperature  box. 


rest  of  view  factors  to  satisfy  the  first  Kirchhoff’s  law,  Eq.  (25): 

Y,f‘j  =  1  (25) 

The  view  factors  are  easily  normalized  to  the  sum  of  those 
considered  in  order  to  satisfy  Eq.  (25).  Eq.  (26)  shows  the  nor¬ 
malized  Fab  for  the  previous  case  of  Fig.  13. 


^AR 


Fab 

Fab  +  Faa  +  2Fad  +  2Fac 


(26) 


When  applying  this  second  method,  care  must  be  taken 
to  evaluate  how  much  radiation  is  being  dismissed  as  the 
errors  made  can  be  not  negligible.  However,  any  of  these  two 
approaches  proposed  is  not  usually  applied  separately  and  the 
most  convenient  method  to  improve  accuracy  without  losing 
robustness  in  the  solution  process  is  to  apply  both  of  them  in  the 
order  exposed  here.  Results  shown  below  have  been  obtained 
this  way. 

A  comparison  between  the  two  models,  simple  and  com¬ 
plex,  described  previously  has  been  made  in  order  to  assess 
their  impact  on  predicted  temperatures.  Bearing  in  mind  that 
the  so  called  complex  model  doubles  the  number  of  equations 
and  unknowns  needed  to  be  solved,  it  is  important  to  know  if  it 
is  worth  increasing  the  complexity  of  the  model,  from  the  point 
of  view  of  improving  the  accuracy  of  results. 

The  analysis  has  been  done  in  two  steps,  noting  that  every 
run  of  the  model  has  been  done  on  the  assumption  of  variable 
convective  heat  transfer  coefficients.  Firstly,  both  models  have 
been  applied  to  the  same  boundary  conditions,  i.e.  air  and  fuel 
flows;  results  are  shown  in  Table  6  in  a  similar  way  to  Table  4. 
Differences  between  both  cases  are  more  important  than  for  the 
convective  analysis  and  range  from  1 .4%  for  the  power  produced 
to  5.5%  for  predicted  mean  and  maximum  temperatures.  Not 
only  this,  the  fact  that  temperatures  calculated  with  the  simple 
model  are  lower  than  those  from  the  more  complex  model  makes 
it,  although  faster,  more  unsafe  to  use  the  former  from  the  point 
of  view  of  mechanical  integrity. 

The  second  step  of  the  analysis  has  been  done  assuming  that 
voltage  and  mean  operating  temperatures  must  be  kept  at  the 
desired  0.65  or  0.35  V  and  900  °C,  respectively.  For  the  latter,  the 
air  flow  entering  the  cell  must  be  reduced  when  the  simple  model 
is  applied.  Results  are  shown  in  Table  7,  where  air  utilization  is 
now  added  in  order  to  evaluate  the  reduced  air  flow  needed  to 
keep  the  operating  temperature  at  the  same  level  when  the  simple 
model  is  used.  It  can  be  seen  that  the  difference  is  not  negligible 
for  low  current  densities  where  fuel  cells  usually  operate. 


Axial  coordinate  (cm) 


Fig.  14.  Predicted  temperatures  with  simple  and  complex  radiation  model.  ASP: 
air  supply  pipe;  PEN:  positive-electrolyte-negative;  0.65  V,  80%  t/f. 


According  to  results  shown  in  Tables  6  and  7,  if  global  per¬ 
formance  of  the  fuel  cell  is  to  be  evaluated,  there  is  no  need  to 
use  the  complex  method  to  describe  radiative  heat  transfer.  The 
infinite  parallel  walls  assumption  works  well  and  errors  are  kept 
below  5%.  However,  evaluating  local  temperatures  can  be  very 
important  for  some  researchers  devoted  to  improve  the  internal 
performance  of  the  cell.  In  this  case,  the  simple  model  cannot  be 
further  used.  Fig.  14  depicts  the  difference  in  predicted  internal 
temperatures  when  using  both  methods  with  the  same  bound¬ 
ary  and  initial  conditions.  Two  aspects  must  be  emphasized. 
Predicted  temperatures  are  higher,  giving  way  to  a  better  perfor¬ 
mance  in  terms  of  power  produced,  and  more  stable,  i.e.  they  do 
not  vary  very  much  from  the  entrance  to  the  exhaust,  when  the 
complex  model  is  used.  This  fact  is  close  to  the  assumption  by 
Haynes  [12],  who  claim  that  the  solid  structure  of  the  cell  be  at 
constant  temperature  along  the  tube. 


6.  Conclusions 

The  work  presented  here  is  based  on  another  work  previously 
published  by  the  authors  [5].  It  focuses  on  heat  transfer  mod¬ 
elling  inside  tubular  fuel  cells  and  tries  to  improve  the  weakest 
aspects  of  those  models  being  used  by  other  authors  currently 
[2,12,13].  The  following  particular  conclusions  with  respect  to 
the  models  presented  can  be  drawn: 
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1 .  If  the  model  is  intended  to  predict  the  global  performance  of 
the  fuel  cell  and  no  internal  information  is  needed,  complex 
models  are  not  necessary.  In  other  words,  a  constant  convec¬ 
tive  heat  transfer  coefficient  with  a  simple  radiation  model 
can  be  used. 

2.  The  consideration  of  constant  or  variable  convective  heat 
transfer  coefficients  has  a  major  impact  on  temperature  at 
both  ends  of  the  cell  tube  but  does  not  affect  the  shape  of  the 
temperature  curve  elsewhere. 

3.  Despite  conclusion  1,  if  internal  temperatures  need  to  be 
evaluated,  the  simple  radiation  model  is  not  acceptable  as  it 
underestimates  temperatures  by  around  40  °C  at  usual  oper¬ 
ating  voltages. 

Finally,  some  remarks  must  be  done  about  the  results  shown 
before.  Fig.  6  shows  the  accuracy  of  correlations  to  evaluate  Nu 
with  respect  to  experimental  data.  For  the  most  accurate  of  them, 
Gnielinski  or  Schleiser-Rouse,  there  is  still  a  10%  uncertainty  in 
the  calculated  value  which  cannot  be  avoided.  Thus,  when  the 
effect  of  using  constant  or  variable  Nu  is  discussed  in  Section  4, 
it  must  be  taken  into  account  that  deviations  of  one  method  with 
respect  to  another  are  added  to  that  uncertainty  of,  at  least  10%. 
In  fact,  bearing  in  mind  that  Re  is  closer  to  its  critical  value,  the 
uncertainty  about  whether  the  flow  is  laminar  or  turbulent  must 
also  be  considered. 

There  is  still  another  effect  that  might  be  thought  to  affect  the 
calculation  of  convective  heat  transfer  coefficient,  the  porosity 
of  the  wall.  This  effect  is  already  considered  when  the  mass  and 
heat  balance  equations  are  applied  to  each  volume  of  a  slice. 
Thus,  the  amount  of  energy  leaving  the  volume  with  the  gas 
which  is  diffusing  through  the  porous  walls  is  included  in  the 
calculations. 


These  two  issues  addressed  here  do  not  affect  significantly 
to  the  discussions  presented  in  previous  sections  of  the  work. 
Comparisons  between  results  from  different  models  just  add 
deviations  from  a  standard  uncertainty  which  is  out  of  the  scope 
of  this  work  and  is,  actually,  very  difficult  to  eliminate. 
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